Structure and dynamics of a model glass: 
influence of long-range forces 
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We vary the amplitude of the long-range Coulomb forces 
within a classical potential describing a model silica glass and 
study the consequences on the structure and dynamics of the 
glass, via molecular dynamics simulations. This model allows 
us to follow the variation of specific features such as the First 
Sharp Diffraction Peak and the Boson Peak in a system going 
continuously from a fragile (no Coulomb forces) to a strong 
(with Coulomb forces) glass. In particular we show that the 
characteristic features of a strong glass (existence of medium 
range order, bell-shaped ring size distribution, sharp Boson 
peak) appear as soon as tetrahedral units are formed. 

PACS numbers: 61.43. Fs - 61.43.Bn - 63.50,+x 



I. INTRODUCTION 

The connection between the structure and the dynam- 
ics of a given system is of particular interest in glasses 
since the structural entities at the origin of some specific 
vibrational features like the Boson Peak (BP) are not 
yet well identified. For instance in silica, the archetype of 
network forming glasses, the existence of the BP has been 
connected alternatively to the presence of SiC>4 tetrahe- 
dra [Jjj , " domains" M or voids || . In molecular dynamics 
(MD) simulations the most important ingredient is the 
force field used to mimic the interparticle interactions: 
the aim is to reproduce as closely as possible the ex- 
perimental vibrational data and then, by looking at the 
atomic positions, to give an interpretation in terms of 
structure. In that sense the so-called "BKS" potential 
H has been widely used to describe silica and the re- 
sults obtained on the structure Pll, the dynamical fij] 
and the thermal properties [|3| compare well with exper- 
iments. This pairwise classical potential is the combina- 
tion of a short-range term and the long-range Coulomb 
interaction (CI) which is obviously the most important 
ingredient of the force field M. One remarkable prop- 
erty of this potential is that CI forces alone are sufficient 
to reproduce the tetrahedral covalent bonding network 
of the silica atoms without the addition of extra terms 
such as angular terms. Our idea in this paper is to vary 
the amplitude of the CI (this is not strictly equivalent 
to varying the range) in order to study the variation of 
the structure and the dynamics of what can be seen as 
a toy silica glass. With this procedure we can study the 
properties of a system going from a strong network glass 



(the silica glass described by the usual BKS potential) 
to a fragile two-component soft sphere glass when the 
CI is turned off. The question we would like to address 
is: how do the structure and the dynamics change along 
this path ? The aim is to study the connection between 
structural features like the First Sharp Diffraction Peak 
(FSDP) or the existence of a tetrahedral network and 
the BP, this connection being the topic of contradictory 
interpretations in the literature fplEoUTJ 



II. SIMULATIONS 

We performed molecular dynamics simulations for mi- 
crocanonical systems containing 216 silicon and 432 
oxygen atoms confined in a cubic box of edge length 
L = 21.48 A, which corresponds to a mass density of 
« 2.18g/cm 3 very close to the experimental value of 
2.2g/cm 3 . Periodic boundary conditions were used to 
limit surface effects. The particles interact via the BKS 
potential whose standard functional form is, for two par- 
ticles i and j (which can be either Si or O): 
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where r%j is the interparticle distance, e the charge of an 
electron and the parameters Aij, Bij and Cij are fixed 
as follows: A Si0 = 18003.7572 and A 00 = 1388.773 
eV;B Sl0 = 4.87318 and B 00 = 2.76 A" 1 ; C Sl o = 
133.5381 and Coo = 175.0 eVA 6 . The partial charges 
used to compute the standard CI are given by qsi = 2.4 
and qo = —1.2. This original form contains an unphysi- 
cal divergence at very short distances which we have cor- 
rected by adding a short range repulsive term (~ 1/rfj). 
In order to vary the strength of the CI (which is cal- 
culated using the Ewald summation method) we simply 
multiply the partial charges by a parameter A which is 
fixed at the beginning of each simulation: 
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A simulation starts from a liquid configuration at 7000K 
which is equilibrated during 50000 time steps (35 ps). 
This system is then cooled to zero temperature at a 
quench rate of 2.3 x 10 14 K/s which is obtained by remov- 
ing the corresponding amount of energy from the total 
energy of the system at each iteration. The aim here is 
not to obtain fine details of the structure of real silica 



which explains the use of a rather fast quenching rate 
compared to the rates that can nowadays be achieved. 
In any case computer quenches are still far away from 
those used in experiments. The quenched glass samples 
are finally relaxed in the microcanonical ensemble dur- 
ing 100000 time steps (70 ps) after which the temperature 
has reached « 2K. We then average the structural infor- 
mations over the last 20000 time steps to avoid transient 
configurations and we diagonalize the dynamical matrix 
corresponding to the final configuration to obtain the vi- 
brational density of states. 

This procedure has been repeated for eleven values of A 
going from A = 0.1toA = l.l. The ideal "A = 0" sys- 
tem would consist of a soft-sphere like arrangement for 
the oxygen atoms (similar to the systems we have stud- 
ied earlier |l2j) together with a random arrangement of 
the silicon atoms in the interstices of the oxygen struc- 
ture. Because of the weak attractive 0-0 interaction 
present in the short-range part of the BKS potential, a 
small Coulomb interaction is necessary in order to sta- 
bilize this system: we found that a value of A around 
0.065 is sufficient which explains why we started our in- 
vestigation at A = 0.1. For this lowest A value we have 
checked on the radial pair distribution functions as well 
as on other structural characteristics that the structure 
was indeed very close to the ideal arrangement described 
above. On the other hand, strong values of A lead to an 
artificial local densification of the samples which are non 
homogeneous (holes start to form in the structure): for 
this reason the maximum value of A is limited to 1.1. Of 
course the standard BKS results are recovered for A = 1 . 



III. RESULTS AND DISCUSSION 

In order to examine the "quality" of the tetrahedral 
network we first study the tetrahedral O-Si-0 angle 9 
which should be ideally equal to 109.47° in a perfect 
tetrahedron as well as the Si-O-Si angle <f> which measures 
the relative position and orientation of two neighboring 
Si04 tetrahedra. To numerically determine the bond an- 
gle O-Si-0 (resp. Si-O-Si), we determine for each Si (resp. 
O) atom the two nearest O (resp. Si) atoms and we cal- 
culate the angle between the two corresponding segments 
Si-0 (resp. O-Si), the result being averaged over all the 
Si (resp. O) atoms. The variation of 9 and <f> as a func- 
tion of A is plotted in figure 1. For A < 0.5 it is obvious 
that the Si04 tetrahedra are not formed and that the 
network is not connected: the system remains analogous 
to a soft sphere glass up to a threshold value of about 0.6. 
This can also be seen in the inset of figure 1 where the Si 
coordination number is plotted versus A. This quantity 
has been calculated by the integration of the Si-0 radial 
pair distribution functions up to a cut-off distance of 2.3 
A which corresponds to the first minimum for all the val- 
ues of A. This coordination number is greater than 5.5 for 



A = 0.1 and is close to 4 for A ~ 0.6 — 0.7. Nevertheless 
for A > 0.5 9, the tetrahedral angle, is almost constant 
and close to the angle in a perfect tetrahedron (indicated 
by the black arrow) while (j>, the angle between the tetra- 
hedra, is continuously increasing with A. This indicates 
that the CI has a direct influence on the medium-range 
order and once the tetrahedra (the building blocks) are 
created, they form a connected network whose charac- 
teristics depend on the strength of the CI. In agreement 
with a previous study M the mean Si-O-Si angle is close 
to 150° for A = 1.0 (standard BKS description) slightly 
larger than the experimental value of cj) (144° [Of indi- 
cated by the hollow arrow). The calculated value closest 
to 144° corresponds to A = 0.7 but this value is not suffi- 
cient to recover other structural features of vitreous silica 
as we will see below, therefore this is not an argument to 
modify the standard BKS charges. 

A second way to analyze the connectivity of the network 
is to monitor the distribution of n-fold rings: an n-fold 
ring is defined as the shortest closed path of alternating 
Si-0 bonds. Therefore an n-fold ring consists of 2n alter- 
nating Si-0 bonds. The ring distribution as a function of 
ring size is plotted in figure 2 for 3 values of A: 0.2, 0.5 
and 1. For A = 0.2 the ring distribution is typical of an 
open structure with a majority of small rings (the esti- 
mated number of large rings should be taken with a pinch 
of salt since finite size effects due to the limited length 
of the simulation box bias the results). For A = 0.5 the 
distribution becomes peaked for rings with a specific size 
(5 and 6) which extend over distances greater than the 
inter-atomic distances: molecules are formed in the sys- 
tem. Finally for A = 1.0 we find the usual distribution 
of rings, peaked at sixfold rings and almost symmetric 
around n = 6 p4| , typical of the silica network. 
A complementary way of analyzing the structure is to 
compute the static structure factor S(q) as a function of 
A. We used the standard way of calculating S(q) JL5| and 
focused our attention on a common feature in glasses: the 
First Sharp Diffraction Peak. This pre-peak is related 
to medium-range structures in the glass and its origin 
has been the subject of vigorous debate over a number 
of years. In figure 3 we represent the amplitude of the 
FSDP as a function of A and in the inset we show the 
shape of S(q) for A = 0.2, 0.5 and 1 in the vicinity of the 
FSDP, together with the experimental curve |Q. With 
decreasing values of A the amplitude of the FSDP de- 
creases and one can say that for A < 0.4 this feature does 
not exist anymore in S(q). Indeed the residual value of 
the amplitude for small values of A is not significant and 
is an artifact of the Fourier transform as can be seen in 
the inset. Together with the results shown in figure 1, 
this indicates that the FSDP appears when the tetrahe- 
dra (already formed) arrange themselves into a connected 
network. 

All this structural information shows that below a cer- 
tain amplitude of the Coulomb interactions the structure 



of our toy silica glass has lost all the medium or extended 
range features and is similar to a soft-sphere glass. From 
this study it appears that the transition value of A be- 
tween a fragile and a network glass is around 0.5. Of 
course the only proper way to describe the covalcnt bond- 
ing network of silica atoms would be to use quantum 
mechanical calculations and therefore no precise physi- 
cal meaning should be attached to this threshold value 
A = 0.5. This value is linked to the specific form of 
the BKS potential which does not contain 3-body terms 
and in which the Coulomb interactions are sufficient to 
build the SiC-4 tetrahedra. Nevertheless the good qual- 
ity of this "simple" potential has been demonstrated in 
a recent study which has shown that SiC>2 glass sam- 
ples generated with the BKS potential are good starting 
configurations for ab initio quantum mechanical molec- 
ular dynamics calculations |15[| . One interesting conse- 
quence of our results is that one should be cautious when 
screened Coulomb interactions are used to simulate real 
silica glasses. Indeed screened interactions are in a sense 
similar to "weak" interactions like the ones we have sim- 
ulated here with the small values of A and mistreat the 
sum in reciprocal space which is made in the Ewald sum- 
mation method. We show here that the screening should 
not be too strong in order to simulate a realistic silica 
glass. 

In parallel with the study of the structure we have 
calculated the vibrational spectrum as a function of A. 
During the diagonalization of the dynamical matrix we 
did not calculate the eigenvectors (for computer time rea- 
sons) but used the eigenvalues (always positive) to com- 
pute the vibrational density. We were of course interested 
in the Boson Peak which appears to be a characteristic 
of glasses even though it can also exist in some crystals 
J2J . This peak, seen in Raman and inelastic neutron scat- 
tering spectra, reflects an enhancement in the density of 
states compared to the Debye distribution and is located 
at frequencies around vb — 1.5 THz in silica. When 
varying A, the energy (and the pressure since we consider 
a microcanonical ensemble) of our system varies appre- 
ciably and this is reflected in the variation of the largest 
frequency of the whole vibrational spectrum, f max , which 
increases from f max « 20 THz for A = 0.1 to ^ max « 40 
THz for A = 1.0. As expected, also the shape of the 
spectrum varies with A: for A small we have basically 
one band since the oxygen and silicon atoms are almost 
identical whereas for large values of A we observe a clear 
separation between an acoustic and an optical band typ- 
ical of the formation of well defined Si-0 bonds. But 
the Boson peak is in all cases clearly present. However 
due to the A dependence of the whole spectrum, it is 
not legitimate to compare the BPs directly and therefore 
we have first determined the normalized density of states 
g(f) using reduced frequencies / = ^/^ ma x- The normal- 



ization is such that J g(f)df = 1. Then we have used 
a standard way to extract the BP from the vibrational 
spectrum which is to plot g(f)/f 2 because in the Debye 
approximation g(f) ~ f 2 at low frequency. This curve 
has then been fitted by a "generalized Lorentzian" [ p7[ : 

1 = /o/ " ■ WTW (3) 

which has been used to extract the position of the BP 
and the corresponding height in a straightforward man- 
ner. In figure 4 we show the g(f)/f 2 data together with 
the fitting curves for three typical values of A and in the 
inset we report the position and the intensity of the BP 
as a function of A. As expected, due to the normaliza- 
tion of the <?(/) curves, the variations of the position and 
the intensity have opposite trends. Apart from a first 
decrease of the intensity between A = 0.1 and A = 0.2 
that we cannot interpret simply (at very small values of 
A, there are some artifacts due to the form of the poten- 
tial as indicated above), one observes a net sharpening 
of the BP when A becomes larger than 0.5. This seems 
to indicate that a strong BP is linked to the existence 
of tetrahedral units or to the presence of an FSDP and 
therefore supports previous interpretations for the origin 
of the BP in silica flit] . Moreover the global variation of 
the BP characteristics with A confirms the recent exper- 
imental observation that the BP increases when going 
from fragile to strong glasses p8| . Of course these re- 
sults, based on a comparison between reduced spectra, 
should be taken with care: in addition to the increase of 
the bandwidth when increasing A there is also a defor- 
mation of the whole spectrum which is difficult to take 
into account simply. Finally, if some significance can be 
given to the apparent decrease of the BP intensity with 
increasing A for A > 0.7, one could conclude that, once 
the tetrahedral network is formed, the increase of the 
medium range order leads to a decrease of the BP inten- 
sity. Of course, such a statement needs to be validated 
by further investigations. 



IV. CONCLUSION 

In conclusion, we have performed classical molecular 
dynamics simulations in a toy silica glass model in which 
we have changed the intensity of the long-range Coulomb 
forces. With this procedure it is possible to analyze the 
structure and the dynamics of a fragile glass (without 
Coulomb interaction) and a strong glass (with Coulomb 
interaction) within the same system. We have shown that 
the formation of a tetrahedral network, the existence of 
a First Sharp Diffraction Peak in the static structure fac- 
tor, a characteristic peaked ring size distribution, as well 
as a relatively sharp Boson peak, are correlated to the 
strength of the Coulomb forces: all these features are ab- 
sent in the fragile glasses and appear once the Coulomb 



forces exceed a potential-dependent threshold. Of course 
our study is limited to a very specific model and other 
ways of creating locally ordered units could have been 
explored. Nevertheless it sheds light on what is really 
happening in a system which can be considered as either 
a fragile or a strong glass. By tuning one single parameter 
(which therefore appears to be the predominant ingredi- 
ent in the simulation of glasses) one is able to connect 
structural features like the FSDP to vibrational specifici- 
ties like the Boson peak and give answers to questions 
previously raised by experimental studies. 

Part of the numerical calculations were done at the 
computer center CINES (Centre Informatique National 
de l'Enscigncment Superieur) in Montpcllicr. 
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FIG. 1. (a) Variation of 9 (O-Si-0 angle : O ) and (f> 
(Si-O-Si angle : •) as a function of A. The filled and hollow 
arrows indicate the perfect tetrahedral value of 9 and the ex- 
perimental value of <f> respectively. In inset the Si coordination 
number is shown as a function of A. 
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FIG. 2. Ring size distribution for A = 0.2, 0.5 and 1. 
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FIG. 3. First Sharp Diffraction Peak amplitude as a func- 
tion of A. In inset the structure factor S(q) is shown in the 
vicinity of the FSDP for A=0.2 (long dashed line), A=0.5 
(dashed line), A=1.0 (thin solid line) and compared to the 
experimental S(q) (bold solid line) for SiCh- 
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FIG. 4. The Boson peak represented as a plot of g(f)/f 2 
versus /, where g(f) is the normalized density of vibrational 
states and / the "reduced frequency", taken equal to i//i/ max , 
where tWx is the maximum frequency of the spectrum. A, 
□ and O correspond to A = 0.2, 0.6 and 1.0 respectively. 
Dot-dashed, dashed and continuous lines are the correspond- 
ing fits using formula (3). In inset the position (O) and the 
intensity (•) of the Boson peak are shown as a function of A. 



